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V  I.  INTRODUCTION 

•  I- 

Just  after  World  War  II  Dr.  R.  Beuwkes  [1]  of  Watertown  Arsenal 
introduced  the  senior  author  tcr^the  theory  of  elastic  stresses  in  thick- 
walled  cylinders  and  the  technology  of  shell-pushing  tests.  The  research 
at  Watertown  Arsenal  culminated  in  the  publication  of  the  Thick  Walled 
Cylinder  Handbook  [2],  a  monument  to  skill  in  analysis  and  computations  on 
a  desk  calculator.  ^Investigation  of  non-metallic  rotating  bands  at  BRL  led 
to  renewed  interest  in  this  area.  We  found  that  we  were  unable  to  inter¬ 
polate  in  the  tables  cited  above  due  to  their  limited  accuracy;  moreover 
the  value  of  Poisson's  ratio  used  in  the  computations  was  not  appropriate 
for  modern  gun  steels.  An  independent  investigation  was  initiated,  using 
residue  theory  in  place  of  Fourier  series  [iy4^5].  Since  the  eigenvalues 
were  complex,  we  required  a  subroutine  for  Bessel  functions  of  integral 
order  and  complex  argument.  The  required  subroutine  was  developed  at  BRL. 

A  Gauss  continued  fraction  was  used  to  reduce  round  off  errors  inherent  in 
series  calculations  [6^7}_ — A  code  giving  accurate  stresses  on  the  outside 
of  the  gun  tube  was  developed  [&y9] .  This  code  was  used  to  calculate 
strains  in  a  highly  instrumented  gun  tube.  The  appropriate  value  of 
Poisson's  ratio  and  Young's  modulus  was  obtained  from  the  Benet  Laboratory. 
Agreement  between  theory  and  measurement  was  good  at  low  velocities,  but 
systematic  deviations  were  observed  at  high  velocities.^  This  result  was 
forecast  in  an  early  paper  by  G.S.  Taylor  [10],  who  used^a  dynamic  version 
of  the  Winkler  theory  for  thin-walled  tubes,  but  his  results  were 
apparently  ignored  by  the  Army.  A  program  based  on  scalarynd  vector  wave 
functions  was  initiated  at  BRL.  The  computations  are  difficult  except  for 
torsion,  which  we  discuss  below.  The  theoretical  work  shows  that  the 
equilibrium  stress  distribution  is  obtained  when  the  velocity  of  travel 
approaches  zero  in  the  limit,  as  one  would  expect  on  physical  grounds. 
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Recently  we  have  resolved  difficulties  in  calculating  stresses  near  a 
discontinuity  of  loading  on  the  inner  surface  of  the  cylinder.  A  method  of 
calculating  elastic  stresses  in  thin-walled  cylinders  was  also  derived,  so 
we  are  able  to  use  the  same  mathematical  formulation  for  wall  ratios 
ranging  from  .01  to  5.  Both  of  these  problems  required  asymptotic  methods 
and  involved  large  values  of  the  complex  transform  variable  in  the  analysis 

II.  NUMERICAL  DIFFICULTIES 


Formulation  of  boundary  value  problems  for  the  infinite  hollow 
cylinder  has  followed  traditional  lines  and  is  not  exceptionally  difficult. 
Real  problems  arise  in  the  numerical  evaluation  of  Fourier  integrals  and 
the  generation  of  Bessel  functions  of  the  second  kind  due  to  the  integer 
arithmetic  of  the  digital  computer  and  its  limited  exponent  range.  Memory 
requirements  and  execution  time  are  relatively  modest  for  the  class  of 
problems  under  consideration.  We  have  considered  four  types  of  error  in 
the  course  of  programming  and  numerical  analysis. 

Round  off  error  is  persistent  and  insidious.  It  is  very  severe  in  the 
evaluation  of  Fourier  integrals  by  quadratures  along  the  real  axis  and  was 
the  principle  reason  why  the  calculus  of  residues  was  used  in  preference. 

It  occurred  in  acute  form  in  calculating  Bessel  functions  of  the  second 
kind.  This  difficulty  motivated  our  development  of  the  subroutine  cited 
above.  In  this  paper  we  discuss  round  off  error  occurring  in  the  manipu¬ 
lation  of  asymptotic  series.  Round  off  error  is  also  a  principle  concern 
in  generating  special  functions  by  recursion  formulas,  where  it  arises  in 
connection  with  stability  criteria. 

A  continued  fraction  obviously  can  be  used  only  for  values  of  the 
variable  and  parameter  for  which  division  by  zero  will  not  occur.  We 
finally  are  able  to  prove  that  division  by  zero  would  not  occur  in  the 
portion  of  the  subroutine  using  Gauss  continued  fractions.  Theorems  of 
Bucholz  [11]  and  Hurwitz  [12]  were  required  in  the  proof  [13].  The 
analysis  is  closely  related  to  Hurwitz  stability  theory. 

Serious  truncation  error  has  occurred  only  in  evaluating  residue 
series  for  the  inner  radius  at  points  very  close  to  the  discontinuity  of 
loading.  Only  recently  have  we  found  a  method  for  improving  the  conver¬ 
gence  of  the  residue  series. 


III.  STRESSES  NEAR  A  DISCONTINUITY  OF  LOADING  FOR  AXIALLY  SYMMETIC  STRESSES 

For  brevity  we  consider  only  axial  stresses  produced  by  a  step 
function  of  pressure  or  shear  applied  to  the  inner  cylindrical  surface. 

The  analysis  of  tangential  stresses  is  similar.  We  superimpose  a  constant 
stress  and  a  discontinuity  stress  to  obtain  the  step  function.  For 
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pressure  loading,  we  have 

t  =0,  0  -  h  on,  r  =  a  and  (la) 

rz  r  u 

Trz  =  °'  °r  =  -  *  V  z<0‘«  ar  =  15  V  Z>0»  r  =  a'  (lb) 

respectively,  and  for  shear  loading 

°r  “  °*  arz  =  *5  V  r  =  a  (2a) 

ctz  =  0,  tr2  =  -  %tq,  z<0;  Tr2  =  %t0,  z>0,  r  =  a  (2b) 

In  both  cases 

ar  =  °’  Trz  =  °«  r  =  b*  (3) 

The  solutions  corresponding  to  (la)  and  (2a)  can  be  obtained  by  elementary 
methods  and  will  not  be  considered  here.  The  discontinuous  stresses  in 
(lb)  and  (2b)  are  represented  by  Cauchy  discontinuous  factors  to 
facilitate  solution  by  separation  of  variables. 

The  stresses  are  derived  from  Love's  stress  function  [14]  in  the  form 
♦  =  [Alg(sr)  +  BKQ(sr)  +  Csrl^(sr)  +  DsrK^ (sr) ] cos (sz)  (4) 


for  pressure  loading  and 

♦  =  [Alg(sr)  +  BKq(st)  +  Csrlj(sr)  +  DsrKj(sr)]sin(sz)  (5) 

for  shear  loading.  If  the  boundary  conditions  are  homogeneous,  we  obtain 
four  homogenous  linear  equations  which  are  satisfied  only  if  the  determi¬ 
nant  of  the  coefficients  is  equal  to  zero.  In  the  case  of  shear  loading, 
we  obtain  from  Eq.  (1)  and  a  number  of  intermediate  calculations  the 
characteristic  equation 

A  =  0 
c 


(6) 
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where 


yp) 

yp) 

y^p) 

«1K1CP) 

yp) 

-yp) 

pi0(p) 

-PK0(P) 

yq) 

yq) 

Bjijfq) 

BjKj (q) 

yq) 

-yq) 

qyq) 

-qyq) 

and 


p  *  sa,  aj  »  [p  ♦  (2  -  2v)/p] ,  q  =  sb,  =  [q  +  (2  -  2v)/q] . 


(7) 


(8) 


The  shear  loading  leads  to  the  characteristic  equation 


As(s)  -  0  (9) 

where 


As(s)  a  -  Ac(s)  (10) 

and  obviously  has  the  same  characteristic  roots. 

The  characteristic  roots  in  the  first  quadrant  of  the  complex  s  plane 
have  the  approximate  value 


s  =  t  / (b-a) 
n  n  v 


(11) 


where 


t_  ~  log  [(2n-l)u]  +  i  (n-h)ff,  n>l 

n  c 


(12) 


The  approximate  values  of  sn  obtained  from  Eq.  (11)  are  improved  by  Newton’s 
method  in  the  complex  plane.  It  should  be  observed  that  all  the  determi¬ 
nants  occurring  in  the  analysis  are  analytic  functions  of  s  even  though 
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logarithms  occur  in  the  series  for  the  modified  Bessel  function  of  the 
second  kind. 

We  find  that 


a 

z 


fo  J  sin(sz)ds/s 

0 


(13) 


is  the  solution  corresponding  to  pressure  loading,  Eq.  (lb),  and 


T 

rz 


jr 


00 


(A2/As)cos(sz)ds/s 


(14) 


is  the  solution  to  the  shear  problem,  Eq.  (2b) .  The  determinants  A^  and  A0 
are  given  by 


0 

0 
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-2K0(p) 

+  pKx(p) 

ix(p)/p 

-K1(p)/p 
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-3K0(p) 
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BjKj (q) 
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where 


a2  =  (2-2v)/p, 


B2  »  (2-2v)/q 


(17) 
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We  obtain  asymptotic  approximations  of  the  integrands  by  using 
Wronskian  relations  connecting  Iq(x),  I^(x),  Kq(x),  and  Kj (x) ,  where  x  =  p 

or  x  *  q,  and  the  leading  terms  of  the  Hankel  asymptotic  expansions  [15]. 
The  leading  terms  are 


!000  * 
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We  find 


(18) 

(19) 


(20) 

(21) 


A1/(sAc)  *  s/(s+sQ), 


A2/(sAs)  =  -  s/(s+sQ) 


(22) 


where 


sQ  *  (7-8v) (b-a)/4ab 


(23) 


On  combining  Eqs.  (13),  (14),  and  (22)  we  find  the  resulting  integrals  can 
be  expressed  in  terms  of  sine  and  cosine  integrals  [16] .  The  approximation 
for  small  z  follows  from  the  fact  that  large  values  of  the  transform 
variable  s  correspond  to  small  values  of  the  argument  z,  according  to  the 
usual  theory  of  Fourier  integrals.  Let  y  be  Euler's  Constant  in  this 
context  and  let  Zq  *  1/Sq  be  a  characteristic  length.  Then,  when  z  is 
positive  and  very  small ,  we  have  approximately 


a 


z 


(24) 


for  pressure  loading  and 
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T 

rz 


IT 


[Y  +  loge(z/z0)] 


(25) 


for  shear  loading.  We  represent  the  logarithm  by  an  integral  of  Fourier 
type  [17,  18].  We  subtract  these  dominant  terms  from  the  integrals  given 
in  Eqs.  [13]  and  [14].  We  obtain 


-  Ac)sin(sz) 
'Ac(s) 


/ 


A2cos(szq) 

sAc(s) 


+ 


CD 


/ 


(4  A2-As) (2cos(sz) -2cos(sz0)ds 

s  As(s)  _ 


(26) 


(27) 


The  integrals  in  Eqs.  (26)  and  (27)  are  more  rapidly  convergent  than  the 
original  integrals  in  Eqs.  (13)  and  (14)  and  will  lead  to  more  rapidly 
convergent  residue  series  when  the  limits  of  integration  are  taken  between 
-  «  and  a".  The  integrals  must  be  re-written  in  exponential  form  as 
illustrated  in  the  torsion  problem  to  insure  convergence  of  the  contour 
integrals. 


VI.  ELASTIC  STRESSES  IN  THIN  WALLED  CYLINDERS 

We  observe  from  Eqs.  (11)  and  (12)  that  the  eigenvalues  of  high  order 
for  a  thin-walled  cylinder  become  very  large  in  absolute  value.  Exponen¬ 
tial  over-run  then  occurs  when  we  use  the  Hankel  asymptotic  series  to 
evaluate  the  various  determinants.  Moreover,  when  we  use  Laplace's 
reduction  of  the  determinant  in  Eq.  (7),  we  find  expressions  like 
p[I»2(p) -i^2(p)  an{j  p[Kg2(p) -K^(p) ]  occur,  together  with  similar  express¬ 
ions  involving  q.  When  these  expressions  are  evaluated  by  means  of  the 
Hankel  asymptotic  expansions,  the  leading  terms  are  cancelled  by  subtrac¬ 
tion,  leading  to  increasingly  severe  round  off  error  as  the  wall  ratio 
approaches  one.  To  overcome  these  difficulties,  we  obtained  asymptotic 
expansions  of  these  expressions  in  which  the  subtraction  occurs  algebrai¬ 
cally  rather  than  numerically.  The  exponentials  were  also  combined 
algebraically,  thus  eliminating  exponential  over  run  for  the  range  of  wall 
ratios  of  interest. 
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Let  [19] 


w  =  AWj  +  BW2  +  Cw^ , 

where 

Wj  =  ir[pl02(p)-pl12(p)l 

w2  “  [PIo(P)K0(p)+pIl(p)Kl(p)i 

w3  =  [pK02(p)-pK12(p)]/ir, 

Then  w  satisfies  the  following  differential  equation. 

3  2  2 

p  w''1  +  2  p  w"  -  (4p  +p)w'  +  w=0 

We  find 

_  ~n 
w  *  I  a  p  , 

‘  0  n 

where  the  odd  numbered  coefficients  are  zero,  a^  =  1,  a2  =  -  j 
an=(n3  -  5n2  +  7n-3)/(4nan_2) 
for  n>2  and  even. 

We  let  w^  »  e2p  w^,  w^  =  e*2p  w_.  Then 

p^Wjf**  +  (6p3  +  2p2)w^"  +  (8p^+8p2-2)w1  ’  *  (8p^-2p+l)  w. 


with  a  similar  equation  for  w,.  We  find 


ELDER,  WALBERT,  ZIJWERMAN 


W,  »  l  b  p‘n  (36) 

i  x  n 

where  b1  =  1,  =  1/8,  and 

(8n-8)  bn  =  (6n2-14n+6)  bnl  -  (n3-Sn2+7n)bn_2,  n>2  (37) 

The  function  w-j  was  treated  in  a  similar  manner.  These  formulas  were 
programmed.  We  obtained  500  eigenvalues  for  a  series  of  wall  ratios 
ranging  from  .01  through  5,  and  the  corresponding  stresses  at  the  outside 
radius,  where  the  residue  series  is  rapidly  convergent. 


V.  STRESSES  DUE  TO  AN  ACCELERATING  LOAD 

We  outline  a  method  of  analysis  based  on  superposition,  an  eigenvalue 
expansion,  interchange  in  the  order  of  integration,  and  the  evaluation  of 
a  complicated  infinite  integral.  Justification  for  the  various  steps  is 
omitted  for  brevity,  but  will  be  presented  elsewhere  in  due  course. 

We  assume  the  outside  cylindrical  surface  is  free  of  stress,  but  the 
inner  boundary  is  subject  to  a  discontinuous  moving  load. 


xr0  =  0,  r  =  b  (39) 

xr9  =  T0  Fo*‘Z,t^ 

where 

F0U,t)  =  h,  z>T(t)  (40a) 

*  -  h,  z<T(t)  (40b) 


and  T(t)  is  the  travel.  We  assume  the  velocity  is  subsonic,  that  is, 
t(t)<C2  where  C2  is  the  velocity  of  the  shear  wave  in  steel.  In  order  to 
use  separation  of  variables  we  assume 
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00 

F0<-'>  -  T  / 

0 


sin[sg(z,t) ]ds 
s 


(41) 


where 


8(2, t)  =  2  -  T(t) 


(42) 


We  assume  a  solution  of  the  form 


Tre  =  F0(z,t)  [R0(r)  -  IQnR(qn,r)] 


(43) 


where 


R0(r)  =  ta2(b4-r4)]/[r‘‘(b4-a4)] 


(44) 


and 


R(qn,r)  =  A(qn)  [I(qnr)K2(qnb)-K2(qnr)Y2(qnb)] 


(45) 


The  eigenfunction  expansion 


Vr) 


X 

1 


R(qn.r) 


(46) 


can  be  obtained  either  by  the  theory  of  residues  or  the  theory  of 
orthogonal  functions.  We  have  used  both  methods  to  determine  the  Fourier 
coefficients  Qn  and  the  results  agree.  The  qn  are  eigenvalues  obtained 
from  the  characteristic  equation  R(q,a)  =  0,  and  are  purely  imaginary  since 
the  problem  is  formulated  in  terms  of  modified  Bessel  functions. 


We  assume 
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FnU,t) 


•if- 


W  (s.i.t) 


where 


Wn(s ,z,t)=(^(s ,t) sin(sz) -Hn(st)cos(sz) 


On  combining  these  remits  and  substituting  the  value  of  t  thus  obtained 
in  the  differential  equation 

3^t  .  .  3t  4t  3^t  .  ,  3“t 

_ li  +  I  _£ _ £0  _ r9_  1_  _ r6 

3  ?  7  *  1  ">  0  v^yJ 

3r“  3  3r“  r  3z^  c*  3t2 


we  obtain  two  ordinary  linear  inhomogeneous  differential  equations  for  a 
Gn(s,t)  and  H^Cs.t) .  We  solve  by  Duhamels  integral  and  evaluate  Wn(s,z,t). 
Duhamel's  integral  will  appear  inside  the  integral  in  Eq.  (48).  We 
interchange  the  order  of  integration.  We  obtain  on  letting  a»t-t^, 


Fn(=.t)  - 


sin[a/s  -q^] sinBsdsdtj 

c7  s  V 


where  S  =  z  -  T(t^)  in  the  above  equation.  We  differentiate  the  inner 
integral  partially  with  respect  to  B,  evaluate  the  resulting  inner  integral 
by  means  of  a  known  formula,*  and  integrate  with  respect  to  B  to  regain  the 


original  function  Fn(z,t).  We  obtain 


t  a 


F„(z.t)  -  H  Q„  c2  pn2  /  f  Jo  pny«2-s2n.v  JSdtj  f51] 

0  0 

2  2 

Where  pn  =  -  qn  ,  and  is  real  and  positive.  Thus  we  have  two  quadratures 
followed  by  a  summation.  In  practice,  the  order  summation  and  quadratures 
should  be  interchanged  to  reduce  round  off  error. 


•Reference  17,  page  472,  paragraph  3.876,  Eq.  (1) 
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We  can  readily  obtain  the  response  to  a  step  function,  then  a  square 
wave  by  superposition  and  translation.  An  additional  convolution  will 
account  for  variable  torque,  which  must  be  obtained  from  the  dynamics  of 
the  shell. 


VI.  DISCUSSION  AND  CONCLUSIONS 

We  have  suggested  methods  of  improving  the  accuracy  of  calculations 
based  on  classical  analysis  without  using  multiple  precision  calculations, 
and  which  are  thus  suitable  for  a  group  in  engineering  or  applied  mechanics. 
We  have  obtained  formulas  for  stationary  loads,  loads  moving  with  constant 
velocity,  and,  in  the  case  of  torsion,  loads  moving  with  arbitrary 
acceleration.  The  method  presented  here  for  solving  the  acceleration 
problem  has  not  been  found  in  the  literature  and  therefore  requires  careful 
justification.  Additional  analysis  and  considerable  programming  are 
required  to  obtain  codes  for  calculating  the  stresses  and  strains.  Only 
then  will  the  results  be  useful  in  interpreting  strains  obtained  with 
instrumented  gun  tubes.  The  work  is  continuing  with  the  time  and  resources 
available . 
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